home *** CD-ROM | disk | FTP | other *** search
/ The 640 MEG Shareware Studio 2 / The 640 Meg Shareware Studio CD-ROM Volume II (Data Express)(1993).ISO / prog / yamp2.zip / VIRTRAM.CPP < prev    next >
Text File  |  1992-01-16  |  15KB  |  622 lines
  1.  
  2. /*
  3. YAMP - Yet Another Matrix Program
  4. Version: 1.2
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/23/92
  7.  
  8. Copyright(c) Mark Von Tress 1992
  9. Portions of this code are (c) 1991 by Allen I. Holub and are used by
  10. permission
  11.  
  12. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  13. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  14. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  15. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  16. FROM USE OF THIS PROGRAM.
  17.  
  18. */
  19.  
  20.  
  21.  
  22. // #include "virt.h"
  23.  
  24. //////////////////////////////////////////// string functions
  25. strtype::strtype( strtype &str ) // copy constructor
  26.   { int len = MAXCHARS;
  27.     len = strlen( str.name );
  28.     len = (len >= MAXCHARS ) ? MAXCHARS-2 : len;
  29.     strncpy( name, str.name, len );
  30.     name[len] = '\0';
  31.   }
  32. strtype::strtype( char *str)
  33.   { int len = MAXCHARS;
  34.     len = strlen(str);
  35.     len = ( len >=  MAXCHARS )? MAXCHARS-2 : len;
  36.     strncpy(name, str, len);
  37.     name[len] = '\0';
  38.   }
  39. strtype::strtype( void )
  40.   {
  41.     name[0] = '\0';
  42.   }
  43.  
  44. strtype strtype::operator+(strtype str)
  45.   {
  46.  
  47.     int len1,len2;
  48.     len1 = strlen(name);
  49.     len2 = strlen(str.name);
  50.     int len = ( len2 + len1 >=  MAXCHARS ) ?  MAXCHARS-len1-3 : len2;
  51.     len = ( len < 0 || len >= MAXCHARS-1 ) ? 0 : len;
  52.     strncat(name, str.name, len);
  53.     name[len+len1]='\0';
  54.     return *this;
  55.   }
  56. strtype strtype::operator+(const char  *str)
  57.   {
  58.     int len1=strlen(name);
  59.     int len2=strlen(str);
  60.     int len = ( len2+len1 >=  MAXCHARS ) ?  MAXCHARS-len1-2 : len2;
  61.     len = ( len < 0 || len >= MAXCHARS - 1 ) ? 0 : len;
  62.     strncat(name, str, len );
  63.     name[len+len1]='\0';
  64.     return *this;
  65.   }
  66. strtype strtype::operator=(strtype str)
  67.   {
  68.     int len=( strlen(str.name) >= MAXCHARS ) ? MAXCHARS-3 : strlen(str.name);
  69.     strncpy( name, str.name, len );
  70.     name[len] = '\0';
  71.     return *this;
  72.   }
  73. strtype strtype::operator=(const char  *str)
  74.   {
  75.     int len=( strlen(str) >= MAXCHARS ) ? MAXCHARS-2 : strlen(str);
  76.     strncpy( name, str, len );
  77.     name[len] = '\0';
  78.     return *this;
  79.   }
  80.  
  81.  
  82. //////////////////////////////////////////////////matrix stack functions
  83. MStack::MStack(void )
  84. {
  85.    static int cnter=0;
  86.    next = NULL;
  87.    stackloc = 0;
  88.    level = 0;
  89.    if ( !cnter ){
  90.       cnter = 1;
  91.       level = 1;
  92.    }
  93. }
  94.  
  95. void MStack::Inclevel( void )
  96. {
  97.    Dispatch->level++;
  98. }
  99.  
  100. void MStack::Declevel( void )
  101. {
  102.  
  103.    (Dispatch->level)--;
  104.    if( Dispatch->level < 0 ){
  105.       printf("ERROR: LEVEL has been decremented too often\n");
  106.       exit(1);
  107.    }
  108. }
  109. void VMatrix::NewReference( VMatrix &ROp )
  110. {
  111.     signature = SIGNATURE;
  112.     r = ROp.r;
  113.     c = ROp.c;
  114.     // release the current header and share it with ROp.hdr
  115.     PurgeVectors();
  116.     v = ROp.v;
  117.     v->nrefs += 1;
  118.     mcheck = v->mm;
  119.   }
  120.  
  121. void MStack::Push( VMatrix &ROp )
  122. {
  123.    ROp.Garbage("Push");
  124.    MStack *temp = new MStack;
  125.    temp->NewReference( ROp );
  126.    temp->Nameit( ROp.Getname(ROp) );
  127.    temp->next = Dispatch->next;
  128.    Dispatch->next = temp;
  129.    temp->level = Dispatch->level;
  130.    temp->stackloc = ++(Dispatch->stackloc);
  131. }
  132.  
  133. VMatrix& MStack::Pop( void )
  134. {
  135.    Garbage("Pop");
  136.    if( Dispatch->next && Dispatch->stackloc){
  137.       MStack *t = Dispatch->next;
  138.       Dispatch->next = Dispatch->next->next;
  139.       delete t;
  140.       (Dispatch->stackloc)--;
  141.       return *Dispatch;
  142.    }
  143.    else {Nrerror( "POP: Stack is empty, cant pop any more\n");}
  144. }
  145.  
  146. void MStack::Cleanstack( void )
  147. {
  148.    while( Dispatch->next->level >=  Dispatch->level
  149.       && Dispatch->next->next    )  Dispatch->Pop();
  150.  
  151. }
  152.  
  153. void MStack::PrintStack( void )
  154. {
  155.    MStack *temp = Dispatch;
  156.    int t= 1 + Dispatch->stackloc;
  157.  
  158.    while( t-- ){
  159.       printf("\n Stack matrix :location level r c:%d %d %d %d: name: ",
  160.       t,temp->level,temp->r,temp->c);
  161.       temp->Showname();
  162.       temp = temp->next;
  163.    }
  164.    printf("\n\n");
  165. }
  166.  
  167. //////////////////////////////////////////////////////matrix class
  168.  
  169. VMatrix::VMatrix( void )
  170. {
  171.    r = c = 1;
  172.    Nameit("t");
  173.    curvecind = 0;
  174.    signature = SIGNATURE;
  175.    SetupVectors(r,c);
  176. }
  177.  
  178. VMatrix::VMatrix( const char *str, int rr, int cc)
  179. {
  180.    r=rr;
  181.    c=cc;
  182.    Nameit( str );
  183.    curvecind = 0;
  184.    signature = SIGNATURE;
  185.    SetupVectors(r,c);
  186. }
  187. VMatrix::VMatrix( VMatrix &ROp ) // copy constructor
  188.   {
  189.           ROp.Garbage("Copy Constructor");
  190.           r=ROp.r;
  191.           c=ROp.c;
  192.           name = ROp.name;
  193.           curvecind = 0;
  194.           signature = SIGNATURE;
  195.           SetupVectors(r,c);
  196.           for (int i=1; i<=r; ++i) {
  197.               for (int j=1; j<=c; ++j) M(i,j) = ROp.m(i,j);
  198.           }
  199.           Dispatch->Cleanstack();
  200.    }
  201.  
  202. VMatrix::~VMatrix( void )
  203. {
  204.    PurgeVectors();
  205.    signature = 0;
  206.    r = 0;
  207.    c = 0;
  208. }
  209.  
  210. //////////////////////////////////// matrix member functions
  211.  
  212. double VMatrix::Nrerror( const char *errormsg )
  213. {
  214.    double x;
  215.    printf("\nMATRIX ERROR: %s \nexiting to system\n", errormsg );
  216.    x = 2;
  217.    exit(1);
  218.    return x;
  219. }
  220.  
  221. void VMatrix::Garbage( const char *errormsg )
  222. {
  223.    int errorint = 0;
  224.    // if( Dispatch->signature != SIGNATURE) errorint++;
  225.    if( signature != SIGNATURE  ) errorint++;
  226.    if( v->mm[-1] != SIGNATURE ){
  227.       printf( "v->mm[-1]= %f\n",v->mm[-1]);
  228.       errorint++;
  229.    }
  230.    if( v->mm != mcheck ) errorint++;
  231.    if( !v->mm ) errorint++;
  232.    if( errorint ){
  233.       Dispatch->PrintStack();
  234.       printf("\nFunction name: %s", errormsg);
  235.       Nrerror("Operating on Garbage");
  236.    }
  237. }
  238.  
  239. void VMatrix::SetupVectors( int rr, int cc)
  240. {
  241.    unsigned int numele = (rr * cc) + 1;
  242.    unsigned int siz =  sizeof(double);
  243.  
  244.    if ( numele > 8190 )
  245.       Nrerror("allocation failure 1, too many doubles");
  246.  
  247.    v = (vdoub *) calloc(1, sizeof( vdoub ));
  248.    v->nrefs = 1;
  249.  
  250.    if ( ! (v->mm =(double *) calloc( numele, siz)) )
  251.       Nrerror("allocation failure 2 in SetupVector()");
  252.  
  253.    v->mm[0] = SIGNATURE;
  254.    v->mm++;
  255.    mcheck = v->mm;
  256.    r = rr;
  257.    c = cc;
  258.  
  259. }
  260.  
  261. void VMatrix::PurgeVectors(void )
  262. {
  263.    Garbage("PurgeVectors");
  264.  
  265.  
  266.    v->nrefs -= 1;
  267.    if( v->nrefs > 0 ) return;
  268.    v->mm--;
  269.    if( v->mm[0] == SIGNATURE ) free( v->mm );
  270.    free( v );
  271. }
  272.  
  273. void VMatrix::DisplayMat(void)
  274. {
  275.    Garbage("DisplayMat");
  276.    int wid = Getwid();
  277.    int dec = Getdec();
  278.    Showname();
  279.    printf( "\n\n" );
  280.    for (int i=1; i <=r; ++i) {
  281.       for (int j=1; j <=c; ++j) {
  282.          printf( "%*.*lf ",wid,dec,m(i,j) );
  283.       }
  284.       printf( "\n" );
  285.    }
  286.    printf("\n");
  287. }
  288. void VMatrix::InfoMat( void )
  289. {
  290.    Garbage("InfoMat");
  291.    printf( "\n" );
  292.    Showname();
  293.    printf( "\nr c: %d %d\n", r, c );
  294. }
  295.  
  296. void VMatrix::Replace( VMatrix &ROp )
  297. {
  298.    ROp.Garbage("Replace");
  299.    if ( r !=ROp.r || c !=ROp.c) {
  300.       PurgeVectors();
  301.       SetupVectors( ROp.r, ROp.c);
  302.    }
  303.    name = ROp.name;
  304.    for (int i=1; i<=r; ++i) {
  305.       for (int j=1; j<=c; ++j) M(i,j) = ROp.m(i,j);
  306.    }
  307. }
  308.  
  309. void VMatrix::operator= (VMatrix &ROp)
  310. {
  311.    Garbage("operator =");
  312.    ROp.Garbage("operator =");
  313.    Replace( ROp );
  314.    Dispatch->Cleanstack();
  315. }
  316.  
  317. ////////////////////// end matrix member functions
  318.  
  319. /////////////////// freind functions of matrix class
  320.  
  321. ///------------- addition
  322.  
  323. VMatrix& operator+ (VMatrix &LOp, VMatrix &ROp)
  324.   {
  325.     LOp.Garbage("operator +");
  326.     ROp.Garbage("operator +");
  327.     if (LOp.r==ROp.r && LOp.c==ROp.c) {
  328.     VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  329.     for (int i=1; i<=LOp.r; ++i) {
  330.         for (int j=1; j<=ROp.c; ++j) {
  331.         temp->M(i,j) = LOp.m(i,j) + ROp.m(i,j);
  332.         }
  333.     }
  334.     temp->name = temp->name+LOp.name+"+"+ROp.name+")";
  335.     Dispatch->Push(*temp);
  336.     delete temp;
  337.     return  Dispatch->ReturnMat();
  338.     }
  339.     else ROp.Nrerror ("Mismatched Matrix sizes in addition function\n");
  340.   }
  341. VMatrix& operator+ (double scalar, VMatrix &ROp)
  342.   {
  343.     ROp.Garbage("operator s+M");
  344.     strtype string;
  345.     char buffer[MAXCHARS]={'\0'};
  346.     VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  347.     for (int i=1; i<=ROp.r; ++i) {
  348.     for (int j=1; j<=ROp.c; ++j) {
  349.         temp->M(i,j) = scalar + ROp.m(i,j);
  350.     }
  351.     }
  352.     gcvt(scalar, NDECS, buffer);
  353.     string = buffer;
  354.     temp->name = temp->name+string+"+"+ROp.name+")";
  355.     Dispatch->Push(*temp);
  356.     delete temp;
  357.     return  Dispatch->ReturnMat();
  358.   }
  359. VMatrix& operator+ ( VMatrix &ROp, double scalar)
  360.   {
  361.     ROp.Garbage("operator M+s");
  362.     strtype string;
  363.     char buffer[MAXCHARS]={'\0'};
  364.  
  365.     VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  366.     for (int i=1; i<=ROp.r; ++i) {
  367.     for (int j=1; j<=ROp.c; ++j) {
  368.         temp->M(i,j) = scalar + ROp.m(i,j);
  369.     }
  370.     }
  371.     gcvt( scalar, NDECS, buffer);
  372.     string = buffer;
  373.     temp->name = temp->name+ROp.name+"+"+string+")";
  374.     Dispatch->Push(*temp);
  375.     delete temp;
  376.     return  Dispatch->ReturnMat();
  377.   }
  378.  
  379. ////-------------subtraction
  380.  
  381. VMatrix& operator- (VMatrix &LOp, VMatrix &ROp)
  382.   {
  383.     LOp.Garbage("operator M-M");
  384.     ROp.Garbage("operator M-M");
  385.     if (LOp.r==ROp.r && LOp.c==ROp.c) {
  386.     VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  387.     for (int i=1; i<=LOp.r; ++i) {
  388.         for (int j=1; j<=LOp.c; ++j) {
  389.         temp->M(i,j) = LOp.m(i,j) - ROp.m(i,j);
  390.         }
  391.     }
  392.     temp->name = temp->name+LOp.name+"-"+ROp.name+")";
  393.     Dispatch->Push(*temp);
  394.     delete temp;
  395.     return  Dispatch->ReturnMat();
  396.     }
  397.     else Dispatch->Nrerror ("Mismatched VMatrix sizes in subtraction function\n");
  398.   }
  399. VMatrix& operator- (double scalar, VMatrix &ROp)
  400.   {
  401.     ROp.Garbage("operator s-M");
  402.     strtype string;
  403.     char buffer[MAXCHARS]={'\0'};
  404.     VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  405.     for (int i=1; i<=ROp.r; ++i) {
  406.     for (int j=1; j<=ROp.c; ++j) {
  407.         temp->M(i,j) = scalar - ROp.m(i,j);
  408.     }
  409.     }
  410.     gcvt(scalar, NDECS, buffer);
  411.     string = buffer;
  412.     temp->name = temp->name+string+"-"+ROp.name+")";
  413.     Dispatch->Push(*temp);
  414.     delete temp;
  415.     return  Dispatch->ReturnMat();
  416.   }
  417. VMatrix& operator- ( VMatrix &ROp, double scalar)
  418.   {
  419.     ROp.Garbage("operator M+s");
  420.     strtype string;
  421.     char buffer[MAXCHARS]={'\0'};
  422.  
  423.     VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  424.     for (int i=1; i<=ROp.r; ++i) {
  425.     for (int j=1; j<=ROp.c; ++j) {
  426.         temp->M(i,j) = scalar - ROp.m(i,j);
  427.     }
  428.     }
  429.     gcvt( scalar, NDECS, buffer);
  430.     string = buffer;
  431.     temp->name = temp->name+ROp.name+"-"+string+")";
  432.     Dispatch->Push(*temp);
  433.     delete temp;
  434.     return  Dispatch->ReturnMat();
  435.   }
  436. //------------- unary minus
  437. VMatrix& operator- ( VMatrix &ROp)
  438.   {
  439.     ROp.Garbage("operator -");
  440.     VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  441.     for (int i=1; i<=ROp.r; ++i) {
  442.     for (int j=1; j<=ROp.c; ++j) {
  443.         temp->M(i,j) = -ROp.m(i,j);
  444.     }
  445.     }
  446.     temp->name = temp->name+"-"+ROp.name+")";
  447.     Dispatch->Push(*temp);
  448.     delete temp;
  449.     return  Dispatch->ReturnMat();
  450.   }
  451.  
  452. //-------------- multiplication
  453.  
  454. VMatrix& operator* (VMatrix &LOp, VMatrix &ROp)
  455.   {
  456.     double sum=0;
  457.     LOp.Garbage("operator M*M");
  458.     ROp.Garbage("operator M*M");
  459.     if (LOp.c == ROp.r ) {
  460.     VMatrix *temp = new VMatrix( "(", LOp.r, ROp.c);
  461.     for (int i=1; i <=LOp.r; ++i) {
  462.         for (int j=1; j <= ROp.c; ++j) {
  463.            sum = 0;
  464.            for ( int k=1 ; k<=LOp.c; k++)
  465.                sum  += LOp.m(i,k)*ROp.m(k,j);
  466.            temp->M(i,j) = sum;
  467.         }
  468.     }
  469.     temp->name = temp->name+LOp.name+"*"+ROp.name+")";
  470.     Dispatch->Push(*temp);
  471.     delete temp;
  472.     return Dispatch->ReturnMat();
  473.     }
  474.     else Dispatch->Nrerror ("Mismatched VMatrix sizes in multiplication function\n");
  475.   }
  476. VMatrix& operator* (double scalar, VMatrix &ROp)
  477.   {
  478.     ROp.Garbage("operator s*M");
  479.     strtype string;
  480.     char buffer[MAXCHARS]={'\0'};
  481.     VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  482.     for (int i=1; i<=ROp.r; ++i) {
  483.     for (int j=1; j<=ROp.c; ++j) {
  484.         temp->M(i,j) = scalar * ROp.m(i,j);
  485.     }
  486.     }
  487.     gcvt(scalar, NDECS, buffer);
  488.     string = buffer;
  489.     temp->name = temp->name+string+"*"+ROp.name+")";
  490.     Dispatch->Push(*temp);
  491.     delete temp;
  492.     return  Dispatch->ReturnMat();
  493.   }
  494. VMatrix& operator* ( VMatrix &ROp, double scalar)
  495.   {
  496.     ROp.Garbage("operator M*s");
  497.     strtype string;
  498.     char buffer[MAXCHARS]={'\0'};
  499.  
  500.     VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  501.     for (int i=1; i<=ROp.r; ++i) {
  502.     for (int j=1; j<=ROp.c; ++j) {
  503.         temp->M(i,j) = scalar * ROp.m(i,j);
  504.     }
  505.     }
  506.     gcvt( scalar, NDECS, buffer);
  507.     string = buffer;
  508.     temp->name = temp->name+ROp.name+"*"+string+")";
  509.     Dispatch->Push(*temp);
  510.     delete temp;
  511.     return  Dispatch->ReturnMat();
  512.   }
  513.  
  514. //-------- elementwise multiplication
  515.  
  516. VMatrix& operator% (VMatrix &LOp, VMatrix &ROp)
  517.   {
  518.     LOp.Garbage("operator M%M");
  519.     ROp.Garbage("operator M%M");
  520.     if (LOp.r==ROp.r && LOp.c==ROp.c) {
  521.     VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  522.     for (int i=1; i<=LOp.r; ++i) {
  523.         for (int j=1; j<=ROp.c; ++j)
  524.         temp->M(i,j) = LOp.m(i,j)*ROp.m(i,j);
  525.     }
  526.     temp->name = temp->name+LOp.name+"%"+ROp.name+")";
  527.     Dispatch->Push(*temp);
  528.     delete temp;
  529.     return  Dispatch->ReturnMat();
  530.     }
  531.     else Dispatch->Nrerror ("Mismatched Matrix sizes in elementwise multiplication\n");
  532.   }
  533.  
  534. //------------- division
  535.  
  536. VMatrix& operator/ (VMatrix &LOp, VMatrix &ROp)
  537.   {
  538.     LOp.Garbage("operator M/M");
  539.     ROp.Garbage("operator M/M");
  540.  
  541.     if (LOp.r==ROp.r && LOp.c==ROp.c) {
  542.     VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  543.     for (int i=1; i<=LOp.r; ++i) {
  544.         for (int j=1; j<=ROp.c; ++j) {
  545.         double d = ROp.m(i,j);
  546.         double L = LOp.m(i,j);
  547.         temp->M(i,j) = ( fabs(d) > ZERO || fabs((d-L)/d) < 1.0E-5 ) ?
  548.                L / d :
  549.                ROp.Nrerror(" zero divisor in elementwise division");
  550.         }
  551.     }
  552.     temp->name = temp->name+LOp.name+"/"+ROp.name+")";
  553.     Dispatch->Push(*temp);
  554.     delete temp;
  555.     return  Dispatch->ReturnMat();
  556.     }
  557.     else Dispatch->Nrerror ("Mismatched Matrix sizes in elementwise division\n");
  558.   }
  559.  
  560. VMatrix& operator/ ( VMatrix &ROp, double scalar)
  561.   {
  562.     ROp.Garbage("operator M/s");
  563.     strtype string;
  564.     char buffer[MAXCHARS]={'\0'};
  565.  
  566.     if( fabs( scalar ) < ZERO )
  567.        ROp.Nrerror(" zero divisor in elementwise division");
  568.  
  569.     VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  570.     for (int i=1; i<=ROp.r; ++i) {
  571.     for (int j=1; j<=ROp.c; ++j) {
  572.         temp->M(i,j) = ROp.m(i,j)/scalar;
  573.     }
  574.     }
  575.     gcvt( scalar, NDECS, buffer);
  576.     string = buffer;
  577.     temp->name = temp->name+ROp.name+"/"+string+")";
  578.     Dispatch->Push(*temp);
  579.     delete temp;
  580.     return  Dispatch->ReturnMat();
  581.   }
  582.  
  583. //////////////////////end of friend functions
  584.  
  585. ////////////////////// matrix functions
  586.  
  587. ////////////// utility functions
  588.  
  589. void VMatrix::Writeb(char *fid, VMatrix &mat)
  590.   /* write a matrix in an binary file */
  591.   {
  592.      int i;
  593.      FILE *file;
  594.      double x;
  595.      char name[MAXCHARS]={'\0'};
  596.  
  597.      file = fopen( fid, "wb" );
  598.      if ( !file ) Dispatch->Nrerror("WRITEB: unable to open file");
  599.  
  600.      mat.Garbage("Writeb");
  601.  
  602.      strncpy( name, mat.StringAddress(), 79);
  603.      i = strlen( name );
  604.      fwrite(&i, sizeof(int), 1, file);
  605.      fputs( name, file);
  606.      i = mat.r;
  607.      fwrite(&i, sizeof(int), 1, file);
  608.      i = mat.c;
  609.      fwrite(&i, sizeof(int), 1, file);
  610.  
  611.     for( i=1; i<=mat.r; i++ ) {
  612.        for( int j=1; j<=mat.c; j++)
  613.          x = mat.m(i,j);
  614.       fwrite( &x, sizeof(double), 1, file);
  615.     }
  616.  
  617.      fclose(file);
  618.      mat.M(1,1);   // reset matrix to base pointer
  619.   } /* writeb */
  620.  
  621.  
  622.